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ABSTRACT 



o\ : 

i Aims. In this paper, we study the propagation and stability of nonlinear sound waves in accretion disks. 
£SJ , Methods. Using the shearing box approximation, we derive the form of these waves using a semi-analytic approach and go on to 
' study their stability. The results are compared to those of numerical simulations performed using finite difference approaches 
I such as employed by ZEUS as well as Godunov methods. 
C" ' ■ Results. When the wave frequency is between Q and 2fi (where Q is the disk orbital angular velocity), it can couple resonantly 
' with a pair of linear inertial waves and thus undergo a parametric instability. Neglecting the disk vertical stratification, we 
derive an expression for the growth rate when the amplitude of the background wave is small. Good agreement is found with the 
results of numerical simulations performed both with finite difference and Godunov codes. During the nonlinear phase of the 
instability, the flow remains well organised if the amplitude of the background wave is small. However, strongly nonlinear waves 
break down into turbulence. In both cases, the background wave is damped and the disk eventually returns to a stationary state. 
Finally, we demonstrate that the instability also develops when density stratification is taken into account and so is robust. 
Conclusions. This destabilisation of freely propagating nonlinear sound waves may be important for understanding the compli- 
cated behaviour of density waves in disks that are unstable through the effects of self-gravity or magnetic fields and is likely 
to affect the propagation of waves that are tidally excited by objects such as a protoplanet or companion perturbing a proto- 
planetary disk. The nonlinear wave solutions described here as well as their stability properties were also found to be useful for 
testing and comparing the performance of different numerical codes. 
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1. Introduction 



Waves and oscillations are expected to be found in accretion disks as they can be excited by a very large variety of 
phen omena. For example, tidal excitation can b e caused by the influence of a binary companion (lArtvmowicz fe Lubow 
1994 ) or by a planet orbiting inside the disk ( Goldreich fe Tremainelll979l Il980 : Nelson et al. 2000|). Spiral density 



waves can be excited as a result of gravitational instability when the disk mass is large enough ( Lin fe Shu 1964: 

I — II 1 | 1| — i ^ 

Too mrelll981t IS hu 1992) . It h as also been suggested that waves are responsible for the puzzling phenomenon of QPOs 

(|Katoll200ll lArras et al. 20061). Recent num erical simulations of MHD turbulence resulting from the magnetorotational 

instability (MRI; iBalbus fc Hawlev 1998 ) have also indicated that den sity waves could be excited by the turbulent 

fluctuations of the flow ( Papaloizou et al. 2004 ; Gardiner fc Stone]|2005f) . 

There have been many studies of the properties of the large variety of waves that can propagate in di sks. Most 

of the e arly work was done using simplifying assumptions and restricting the analysis to the linear case (|Lin et aF 



1990allb1 ; iLubow fc Pringlejll993l: iKorvcanskv fc Pringlelll995t iLubow fc Ogilvielll998t lOgilvie fc Lubowlll999fh Some 



other work extended these results to the nonlinear regime, either semi-analytically, but retaining terms up to second 
order in an expansion in powers of the amplitude (jLarson 1990), or numerically in order to study wave dissipation 



( Bate et al.ll2002t iTorkelsson et al.l fcoOO) 



Despite all these efforts, a complete semi-analytic solution to the wave equation, valid in the nonlinear regime for 
large amplitude, is still lacking. This would not only be useful for studying the properties of waves, but would also 
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provide a valuable test for numerical codes that are more and more extensively used in the accretion disk community. 
The goal of this paper is to derive such a solution and to study its properties and stability. Of course, the need to 
derive a solution that can be found using semi-analytic methods imposes strong constr aints: we use an isothermal 
equation of state, adopt the shearing box approximation (jGoldreich fc Lvnden-Belllll965j ) and restrict the analysis to 
the case where there is no variation in the direction of rotation or axisymmctry. But all these approximations simplify 
the analysis which can then be done including nonlinearities. 

In this paper we are concerned with parametric instability occurring for non linear density waves treated in a local 
shearing box approximation. However, parametric instability appears also t o be generic for global non axisymmetric. 



shearing box approximation. However, parametric instability appears also t o be generic lor global n 
non m agnetic disk flo ws such as may occur w hen the disk is either eccentric ( Goodma"nlll993t lOgilviel 



20011 : iPapaloizou 



2005bllaK or warped (jTorkelsson et al. I l2000h . It is als o related to the elliptical instability observed in rotating flows 
that is apparently connected to turbulent breakdown (jKerswellll2002l ). Such turbulent breakdown does seem to occur 
for some of the large amplitude non circular flows studied here. All of this makes the phenomenon an important one 
for study in the general context of the dynamics of accretion disks. Being a local instability that occurs through the 
interaction of a dense spectrum of inertial waves with a primary flow, modes with arbitrarily small length scales, 
depending ultimately on the dissipative properties of the system may be excited, providing a challenge for numerical 
simulations. 

The plan of the paper is as follows. In section [21 we derive the equations describing the propagation of nonlinear 
density waves in the shearing box approximation. The solutions to these equations are compared with the results of 
ID numerical solutions. In section[3j we study their stability. When the frequency of the wave lies in the range [f2, 20], 
where Q is the orbital angular velocity, it can interact resonantly with two linear inertial waves (whose individual 
frequencies, being bounded by 0, comb i ne to equal the wa ve frequency) and be unstable to a parametric instability 
(| Goodman! 1 1993t iTorkelsson et al. 1 120001 |Papaloizoull2005ah . Neglecting vertical density stratification in this section, 
we derive an expression for the growth rate of the instability and describe the properties of the eigenmodes. These 
results are compared in section [4] with 2D numerical simulations of vertically unstratified shearing boxes. Both finite 
difference schemes as used in ZEUS and NIRVANA (see below) as well as Godunov methods are shown to give results 
that agree very well with each other and with the results of the linear analysis. Finally, in section [5j we investigate the 
effect of vertical stratification, showing the occurrence of parametric instability in that case and relating the results 
to the unstratified case. Finally we discuss our results and summarize conclusions in section [6l 



2. Free waves 

2.1. Basic equations 

In common with many studi es of the local dynamics of accr etion disks orbiting about a central mass, we adopt the 
shearing box approximation (jGoldreich fc Lvnden-Belllll965( ). In this approximation, a reference frame rotating with 
a fixed angular velocity ft, which corresponds to the Keplerian angular velocity at a specified radius Rq, is adopted. 
The origin, at the centre of the box, is located at some specified point on the circle of radius Rq and a system of 
Cartesian coordinates (x, y, z) with associated unit vectors (i,j,k) is used. Gravity is due to the central object only 
and an expansion of the gravitational acceleration correct to first order in x and z is adopted. Then the equations of 
continuity and momentum conservation can be written in the form 

dp 
dt 



V-{ P v) = 0, 



i:)v 
~dt 



+ (v-V)v 



2Uxv = — VP 
P 



3tt 2 xi - il 2 zk . 



(1) 
(2) 



where p is the density, v the velocity and P is the pressure which is related to the density by an equation of state 
(EOS). In this paper we restrict consideration to an isothermal EOS, for which 

2 



P 



pc 



(3) 



where c = QH is the isothermal sound speed with H being the scale height. In the absence of waves, the equilibrium 
solution takes the form 

( z 2 

p = p exp 



2iJ 2 ) 



3 



(4) 
(5) 



where po, being the midplane density, is a constant taken to be unity in the following. In some cases we neglect the 
vertical component of the gravitational acceleration, — Sl 2 zk, in Eq. Then vertical stratification is also neglected 
and p = po throughout. 
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2.2. The wave equation 

Initially, for simplicity, we neglect vertical stratification and assume that the state variables do not vary along the 
direction of rotation or the vertical direction and thus depend only on x and t. To derive the equations that describe 
wave propagation in this system, we adopt a Lagrangian description. We define X(Xo,t) to be the position, or x 
coordinate, of a fluid element that would be at the fixed location x — X in the absence of any disturbance. The y 
component of the momentum conservation equation can be written 

Dv v BX 

where the convective derivative is 
D d 8 

m = di + v *d-x> (7) 



and we have used the fact that 
DX 



(8) 



Dt 

Eq. © can be integrated to give 

Vy + 2nx = f(X ) ■ (9) 

where f(Xo) is a function of Xq to be determined. 

Applying the same procedure to the equation of motion in the x direction, we find 

^ + Q 2 x = _ldP +2QfiXoh (10) 
Dt p ox 

where we have used Eq. @ to eliminate v y . In the absence of disturbance, v x = and the pressure is uniform. This 
means we must have 

f(x ) = ~nx . (ii) 

Thus Eq. ([TO]) becomes 

^+ "7* • < 12 > 

In order to express the density in terms of X, we use the conservation of mass. This gives 
dX ( dt \ 

p mb='{ 1 + ak) =po > (13) 

where we have introduced the displacement £(Xq, t) through the relation X(X , t) = X + £(Xq, t). 



2.2.1. Travelling waves 

Since the natural boundary conditions for the shearing box are periodic in shearing coordinates, we consider travelling 
wave solutions that depend only on the phase $ = Xo — Ut, where U is a constant the phase velocity. Then we have 
£ = Using Eq. (JT3J) to eliminate the density, we obtain 

U 2 ?-L + n 2 £- - — (U) 

which is the governing equation that describes travelling wave propagation. It has been derived assuming no magnetic 
field is present. We comment here that a vertical or an azimuthal field (or both) that depend only on <!> can easily 
be incorporated into the above analysis. However, these details are beyond t he scope of this paper. We comment also 
that there are some similarities of the above analysis to that of [Larson He considered travelling waves in a disc 



though not explicitly in a shearing box. However, his analysis focused on spiral waves and only terms of second order 
or lower in the wave amplitude were retained. This is inexact and inapplicable to the large amplitude waves considered 
here. 
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Fig. 1. Profile of the function V(p) for different values of (c/U) 2 . The solid, dotted, dashed, dotted-dashed and 
dotted-dotted-dotted-dashed curves respectively correspond to (c/U) 2 — 0.005,0.02,0.1,0.25,0.5. 



2.3. First integral of motion 

By considering a first integral of Eq. (fT4")) , it is possible to interpret the problem as corresponding to the motion of a 
particle in a potential well. To see this, let p be defined through 

Using this, Eq. (fT4"l) can be written 
- 2 N dp 

(i+p) 2 ;^ 



^-7^2)^ + ^ = 0. (16) 



which can be integrated with respect to £ to give 



where the dimensionless function V(p) is given by 

V(p) = 4 + ic/U? (-^--Hl+p)) (18) 



1+p 

and Uo is a constant. Equation (| 1 T[) shows by analogy with classical mechanics that the evolution of p as a function of 
$ can be periodic, consisting of oscillations in the wells of the function V(p), should they exist. The bounds of such 
oscillations are determined by the condition V(p) = (U0/U) 2 . By definition, £ and p should oscillate around p = 
which should be a minimum of V and be such that p > — 1. Such a situation is only possible if c/U < 1. The form 
of V(p) is shown in figure [T] The different curves correspond to (c/U) 2 = 0.005 (solid line), 0.02 (dashed line), 0.1 
(dotted line), 0.25 (dotted-dashed line) and 0.5 (dotted-dotted-dotted-dashed line). The form of V(p) is similar for all 
these values of (c/U) 2 . It increases for — 1 < p < — 1 + c/U, decreases for —1 + c/U < p < and increases again 
for p > 0. The well allows p to oscillates around 0. Figure Q] also shows that the amplitude of the oscillations of p is 
bounded such that \p\ < 1 — c/U 



2.4. The linear limit 

When \p\ is small, the wave becomes linear. This happens when \p\ -C 1 — c/U. In this limit the potential is parabolic, 
being given by 

V(p)~^-(1-(c/U) 2 ) . (19) 
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Fig. 2. Left panel: radial profile of the density for a wave of velocity (c/U) 2 — 0.3, for different wave amplitudes: 
pmax/po = 1-1 (dashed line), 1.4 (dotted line) and 1.8 (solid line). The solution shows a cusp around its maximum as 
the amplitude of the wave is increased. Right panel: same as the left panel, but for (c/U) 2 = 0.9 and p m ax/po = 1.01 
(dashed line), 1.03 (dotted line) and 1.052. Note the smaller wave amplitude and the smaller periods in that case 
compared to the case (c/U) 2 — 0.3. 



The system then behaves as an harmonic oscillator. In this linear limit, Eq. (|14|) becomes 

(u 2 -c 2 )^ + n 2 z = o. (20) 

The solution is 

£ = Co cos(kX a - ut) + & sin(kX - cot) , (21) 

with k = £l/^/U 2 — c 2 , and u = kU which together imply the standard dispersion relation for linear sound waves 
modified by rotation in the form ui 2 = k 2 c 2 + f2 2 . 

2.5. The effect of vertical stratification 

The same equations as above apply when the disk is vertically stratified. The nonlinear wave has exactly the same 
properties in that case. Both the perturbed and unperturbed density are simply multiplied by a factor exp ^— -^jj^J 
which contains the entire vertical dependence. 

2.6. Solutions of the wave equation 

Using Eq. (|15p . the wave equation (TT4"]) can be written as a system of two coupled first order differential equations: 



d<S> 

dp _ fny (c/u) 2 £, 



p. (22) 

(23) 



d<S> V c J l _ ( c ^ 2 ' 

v 7 1 (i+p) 2 

Once the phase $ and displacement £ are normalized by the natural shearing box length scale Q/c, a solution of the 
wave equation is specified by the parameter c/U and the maximum value of the density in units of po, p max /po- The 
latter fixes the initial conditions for the integration, which is performed until the solution reaches the next maximum. 
This procedure gives the spatial period T w of the wave, or, equivalently, its wavenumber k x , w — 2tt/T w . 

The solutions we obtained are illustrated in figure^ The left panel shows the structure of the wave when (c/U) 2 = 
0.3 for different values of the relative wave amplitude p max / po = 1.1, 1.4 and 1.8 (we note that the maximum possible 
free wave amplitude is Pmax/Po = U/c = 1.825 for these parameters). For these three values, we respectively found 
T w — 9.587?, 9A3H and 9.23-ff. As the amplitude of the wave is increased, nonlinearity becomes more and more 
important and a cusp eventually forms around the location of the density maximum. The right panel shows the 
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Fig. 3. Results of the numerical tests illustrating the propagation of ID nonlinear waves in the shearing box. The 
parameters are such that (c/U) 2 = 0.3 and p max /po — 1-05 (upper row) and p max /po = 1.8 (lower row). The resolution 
is with 320 grid points in the radial direction. Upper and lower left panels show results obtained with ZEUS using a 
Courant number Co = 0.5, the upper and lower central panels are obtained with ZEUS using Co = 0.1 The upper and 
lower right panels are obtained with RAMSES using Co = 0.7. For the upper panels, the plots are for times (measured 
in cycles or wave periods, with one cycle or period being the time for the wave to cross the box) t = (solid line), 
t = 250 (dotted line) and t = 500 (dashed line). For the lower panels, the plots are for times t = (solid line), t = 30 
(dotted line) and t = 60 (dashed line). Note that the profiles at t — are computed by solving the wave equation as 
described in section |2~21 



structure of the wave when (c/U) 2 = 0.9. Again, three different wave amplitudes illustrate the transition from linear 
to nonlinear wave profiles. These are given by p m ax/ Po — 1.01,1.03 and 1.052. We note that similar waveforms occur 
for both values of c/U. However, for corresponding waveforms, when (c/U) 2 = 0.9, the amplitudes and periods of the 
wave are much smaller than when (c/U) 2 = 0.3. 



2.7. Numerical tests 



The results presented in section[2]can be used to provide ID te sts for hydrodynamic co des. Here we pres ent such test cal- 
culat ions performed with two finite differ ence codes, ZE US (jHawlev fc Stondll995h and NIRVANA (jZiegler fe Yorkd 



19971 ). and a Godunov code, RAMSES 
Godunov code is de scribed in detail in 
(jFromang et al.ll2006l ). 



Tevssier l2002h. Th e imp lementation of shearing box simulations using a 
Gardiner fc Stond (|2005h and was subsequently applied using RAMSES 



We consider only the case (c/U) 2 — 0.3 here since this is the case we will study in more detail in section [U We 
perform different ID numerical simulations of wave propagation for different wave amplitudes. In each case, the size 
of the box is set equal to the wavelength of the wave. We used periodic boundary conditions in shearing coordinates 
and followed the wave over many cycles. 

The point of these tests is to demonstrate that we can successfully follow free propagating non linear waves. In 
particular we need to verify that this can be done under the conditions adopted in section [¥] where we study the 
hydrodynamic stability of the waves. For that reason, we present results obtained with a single resolution having 320 
grid cells in the radial direction. In these simulations, the timcstcp At is defined by 



At = C - 



Ax 



max \v. r 



(24) 



where Aa; stands for a grid cell siz e and Cp is the Co urant number, which should be positive and smaller than one for 
the CFL condition to be satisfied ( Press et al. 1986f) . 
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2.7.1. Results 

We considered waves with two amplitudes. These are given by p m ax/po — 1-05 and p m ax/po = 1-8. The results obtained 
with ZEUS and RAMSES are illustrated in figure [3l Results obtained using NIRVANA are very similar to the former. 
Each of the panels in figure [3] shows the radial profile of the density at different times measured in wave periods. The 
upper panels illustrate results obtained for the low amplitude case. The lower panels illustrate results obtained for the 
high amplitude case. In the former case, these are shown for t = 0, t = 250 and t = 500. In the latter case, results 
are shown for t = 0, t = 30 and t = 60. The left panels show results obtained with ZEUS with Co = 0.5, the central 
panels show results obtained with ZEUS using Cq — 0.1 and the right panels show results obtained with RAMSES 
with Cq = 0.7. 



2.7.2. Effect of Courant number 

Firstly, it is apparent that the ZEUS results obtained with Cq — 0.5 are never acceptable for long simulations. When 
the wave amplitude is small, we observed that its profile is dispersed on timescales of hundreds of periods. When the 
wave amplitude is large, small scale oscillations indicative of numerical instability, appear on timescales of tens of 
periods. 

As indicated in the central panels, both problems are cured by decreasing the Courant number. The wave drift 
relative to the centre of the box seen at large times in both cases is due to the modification of the wave period as its 
amplitude is slowly damped by numerical dissipation. 

Finally, as shown in the right hand side panels, the results obtained with RAMSES using Co = 0.7 are very stable 
and less dissipative than the ZEUS results. In conclusion, despite differences in the codes, we found that both Godunov 
and finite difference codes are able to follow the evolution of these waves for a long enough timescale to perform the 
simulations we will present in section |4] provided the Courant number is chosen appropriately. For ZEUS or NIRVANA 
this means that Cq < 0.1. 



3. Stability of free waves 

Having found exact nonlinear solutions corresponding to freely propagating density waves, an important issue is their 
stability. A possibility is that such waves undergo parame tric instability through interaction with lower frequency 
inertial or gravity waves ( Goodma"nl 19931 : Papaloizou 2005al ). In addition to the production of these additional waves, 
such instabilities may result in the production of turbulence and the decay of the initial wave amplitude and so affect 
its ability to propagate over large distances. 

Accordingly we now study the stability of non linear density waves to parametric decay into inertial waves (gravity 
waves being absent in our model). This can be done while maintaining axisymmetry, ie. we assume dependence of the 
state variables on x and z only. 



3.1. Stability analysis 

To study the stability of the nonlinear waves described above we start from the equations of motion in Lagrangian 
form. Following the discussion of section [2] when the vertical component of the gravitational acceleration, and thus 
vertical stratification is neglected, these are seen to take the form 



D 2 X 



Dt 2 



D 2 Z 



Dt 2 



1 dP 
p dx 
IdP 
p dz 



(25) 
(26) 



We now want to study the evolution of a small perturbation of the state defined by the unperturbed wave. To do so, 
we write 



X 
Z 



Zq + & , 



(27) 
(28) 



where £ = (£ x >£z) has a small magnitude in comparison to that of £ x = £ the displacement associated with the wave 
used above. We therefore ha ve AX = £ x and AZ = whe re A stands for the Lagrangian perturbation. Upon noting 



that A and D/Dt commute (|Lvnden-Bell fc Ostrikerlll967l ). we obtain 



D%x 
Dt 2 



o dx 



, ox 



(29) 
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Dt 2 



-A 



1 dP 

p dz 



-A 



dh 
~d~z 



(30) 



where h = c 2 In p is the enthalpy. Writin g the Eulerian variation of a quantity / by /', we have the relation A/ 
/' + £.V/ (|Lvnden-Bell fc Ostrikerlll967l ). from which we obtain 



dz 



dh 1 | - <9 2 fo 



Using the continuity equation /?' + V(p£) = 0, we write 



ti 



c 2 dp^ c 2 dp| 2 



p dx p dz 



Combining Eq. (|29]), Eq. (30]) and Eq. (33]), we get 



d d , 



2 d% 



dx 2 



n2 d2 jz c 2 dp dj x 
dxdz p dx dx 



2 d 2 £ z 2 d 2 £ x c 2 dp d£ x 
c — + c 



dz 2 



dxdz p dx dz 



(31) 
(32) 

(33) 

(34) 
(35) 



where we used Eq. (7]) that defines the convective derivative. Adopting a reference frame that moves with the speed 
U of the background wave (discussed in section l2~2l 12.6ft . we also have p = p /(l +p) and v x = —U{l+p) where 



P = d^ x /dx. These, along with the relation v y 



\ilX — tt^Ci fully determine the system. Finally, note that 



dissipation is neglected in our analysis. In Appendix [5] we describe briefly how to incorporate it and the changes it 
produces. However, in the following, we assume the evolution is inviscid as it simplifies the analysis. 

3.2. The case p = 

In the absence of a background wave, p = and the previous system reduces to the simpler form 

2 Q2f a2<T 



--U— 

dt dx 



l x + n% = c 2 ^ ■ - 2 



1-4 1 <■ 



dx 2 



+ C 



d% 
dxdz 

d% 



(36) 
(37) 



dz 2 dxdz 

In this case, the radial displacement of a fluid element can be taken to have the form 

i=f = (£",£") cx e K->i.nt+nk m , v <'-k.») ( 38 ) 
where n, as required by the radial periodicity condition, is an integer. Then Eq. (|36|) and Eq. (|37p lead to the conditions 



[(wj. 



£1 ti h x w c ]£ x -\- Tik x ^ w k z c ^ z 







(39) 
(40) 



nk x , w k z c 2 & + [(uj, n - ncu) 2 ~ k 2 z c 2 ]C = 
From these a dispersion relation can be derived in the form 

[w/,„ - H 4 ~ ((n 2 k 2 x , w + k 2 z )c 2 + ft 2 )K„ - nuj} 2 + n 2 k 2 c 2 = . (41) 
The low frequency branch of this dispersion relation corresponds to inertial waves for which the frequency is given by 

+ fi 2 ' 



(wr,„ - nui) =uj„ = 



1-4/1- 



Ak 2 n 2 c 2 



(k 2 c 2 + n 2 ) 2 _ 

mber limit, th 

these inertial waves is thus smaller than £1 in the original non-comoving frame. 



(42) 



where k 2 — n 2 k 2 w + fc 2 . In the large wavenumber limit, this reduces to W/ )TS — nuj — ±uj n — • The frequency of 
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3.3. The case p ^ 

When p 0, but is in the small wave amplitude limit (i.e. \p\ <C 1), terms of higher than first order in p may be 
neglected, so that the system of equations becomes 



O H~ 



dx 2 
<9z 2 



+ c 



dxdz 
dxdz 



i dp dtix 
dx dx 

• dpd^x 
dx dz 



(43) 
(44) 



In this limit, the background wave is linear and p can be written p = ecos(fc XjU ,x), with e being a small parameter. 
If the frequency of the background wave, uj, lies between fl and 2f2, then it is possible for two inertial waves with 
wavenumbers nk XiW and (n + l)k x w to interact with each other by coupling through the original wave. Working in 
the frame comoving with the wave, noting that the time dependence is separable, we write the perturbed eigenvector 
as a linear combination of inertial waves in the form 



£ _ fg" _|_ e rj") gii^i ,n+Sa n )t-Hnk X:W x-ik z z _|_ _|_ e ^j" + 1 ) 



,i(wi, „+i+5o-„ + i)t+i(n+l)fe a - 



,jX— tk z z 



(45) 



where the wave frequencies are respectively shifted by Sa n and 8a n+ i which are both taken to be of order e. For the 
time dependence to be separable, so allowing the modes to interact, we require that 



Vi,n + fan = ^Ln+l + 8a n+ i . 

Using Eq. (|42[) . if ojj n = uj n + nuj and Wj, n +i = —u> n+ i + (n + 1)uj, this is equivalent to 
uj n + u) n +i = uj + 5a n +i — 6a n = uj + Suj . 



(46) 



(47) 



Here uj n , uJ n -\-i, uj and thus Suj are given, while (|47p gives an interaction condition involving 5a n and 8a n +i. A further 
condition on the latter two, enabling the completion of the solution, is obtained as follows: 

Four linear equations relating the total of eight components of £ , £ ,fj n and fj can be obtained after substituting 



into Eq. (|4"3"|) and Eq. (|44[) and requiring that the resulting Fourier coefficients corresponding to the terms in Eq. (|4"5 
vanish. Setting Dj = w 2 — £1 2 — j 2 k 2 w c 2 these are 



e [D n fj% + nk XiW k z c 2 ^] + 2u n 5<r n t£ + e(n + 1) 
e [D n+1 ^ +1 + (n + l)k x , w k z c 2 r%+ 1 } - 2uj n+1 5a n+1 ^ +1 - 
e [{w 2 n - k 2 z c 2 )J] n z + nk x<w k z c 2 f%] + 2uj n 8cx n Q + e(n + 1) uJUJ n +i 
[(«&f l ~ k 2 z c 2 )v: +1 + + l)fc^fc z c 2 C +1 ] - 2uj n+l Sa n+1 C z +1 - en 



ww n +i — — 



, ,2 k 2 r 2 
UJ «/_ ,„C 



= 



UJ 

T 



Sz 2 z x i w ^ x 



f n -i- —k h r 2 P n 
Sz i 2 z x l w ^ x 



(48) 
(49) 

(50) 
(51) 



Eq. (|39[) and Eq. (|40[) and their counterparts for n. — ► n + 1 give an additional four linear equations. Being linear, in 
order that the total of eight equations be soluble, their determinant must vanish. This leads to the requirement that 



5a n Sa n+ i 

where 
A n 



e 2 (A n + B n )(A n+1 + B n+1 ) 1 



C n C n 



+i 



(52) 



= [(n + l)(k 2 z c 2 - uj 2 n ) + n(n 2 + (n+ ifk^c 2 uj 2 )] uj (uj n - f ) 
= [n(k 2 z c 2 oj 2 n+1 ) + (n + l)(n 2 + n 2 kl w c 2 OJ 2 )] uj (uj n - |) 



B n +1 

C n +i 



{k 2 z c 2 -{n+l)uj 2 n ) 

2 

-{k 2 z c 2 + nul +1 ) 



n 2 + k 2 c 2 + n 2 k 2 x w c 2 - 2ljI 



Q 2 + k 2 c 2 + (n + l) 2 kl w c 2 -2ujf l+1 
Upon noting that Sa n+ i = Sa n + Suj, we find 



8a 



± 



(53) 



10 



S.Fromang & J.Papaloizou: Nonlinear sound waves propagation 




Fig. 4. Properties of the unstable modes of the parametric instability for a wave characterized by (c/U) 2 = 0.3 and an 
amplitude such that the maximum density p ma x/po ~ 1-05. The two panels respectively show the vertical wavelength, 
in units of the scale height, of the unstable mode and the growth rate in units of the wave period as a function of the 
radial wavenumber n. In both panels, the solid lines represent the large n limit given by Eq. (|55[) and Eq. (|56p . Note 
that the mode vertical wavelengths become dense at large n such that the requirement of being equal to the length of 
a prescribed vertical domain such as the scale height when multiplied by some integer can be approached arbitrarily 
closely. 



When Slu < 2a, there is a non zero imaginary part to both Sa n and 5a n +i and the amplitude of the two coupled 
inertial modes grows: this is a parametric instability. The maximum growth rate a is obtained when 5u = 0, which 
occurs when the resonant condition satisfied by the inertial waves is: 



'■>n+l 



(54) 



When Su) does not vanish, the growth rate of the instability is reduced, as also noted bv lGoodma n (1993). We note that 
on account of symmetry with respect to reflection in the midplane, the growth rate does not change when k z — > —k z . 
This enables solutions to be classified according to their symmetry with respect to reflection in the midplane. 



3.4. Unstable mode properties 

In this section, we study the properties of modes that undergo parametric instability for the special case (c/U) 2 = 0.3. 
We focus here on modes that occur when the resonant condition given by Eq. ([54]) is satisfied exactly. 

For a given background wave amplitude, the spatial period of the background wave is calculated as described 
in section \2.G\ The radial wavenumber of the background wave k XfW and the wave frequency to are then fixed. The 
instability couples two inertial waves with radial wavenumbers nk x ^ w and (n+l)k XtW . These should satisfy the resonant 
condition, given by Eq. (|54|) . and the dispersion relations, given by Eq. (|41[) and its counterpart with n — ► n + 1. For 
given values of n, u) and k x ^ w ; lu u , uj n +\ and k z can be determined from these algebraic conditions. At first we shall 
ignore any additional constraints on k z that may arise from boundary conditions such as the requirement that the 
vertical wavelength should be the scale height divided by an integer. 

We note that the large wavenumber limit (n> 1) the solution is given by 



ui n = uj n+1 = — , k z = nk x<w d 4 ^ 2 _ ji ■ { r,r> ' ! 

In the same limit, Eq. ([53")) for the growth rate gives 



e _ / c \ / w 



2 „ „s A 2 m 2 - u? 



° ^ 8 n [u) In J {k ^ H) V + OW^^ 2 - ) ' (56) 
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Fig. 5. Snapshot of the nondimensional vertical velocity in the (r, z) plane showing the structure of the unstable 
mode of the parametric instability when n = 4. The parameters of the background wave are (c/U) 2 = 0.3 and 

Pmax/PO = 1-05. 



To illustrate the properties of the unstable modes, in figure 0] we plot the growth rate a — <j(2it/uj) and vertical 
wavelength as a function of n for a wave whose amplitude is such that p max /po — 1.05. It is seen that results obtained 
assuming the large n limit applies are always in good agreement with the actual results. We comment that in the limit 
of large n the modes become dense in the sense that the requirement that the product of the vertical wave length 
and some integer be the extent of a vertical domain or the scale height can be satisfied with arbitrary accuracy. This 
follows from Eq. (|55|) and the fact that the rationals are dense. 

Finally, the structure of the eigenmode is illustrated in a particular case in figure [5] This shows a snapshot of the 
vertical velocity in the (r, z) plane for a mode with n = 4, in the case of a background wave for which (c/U) 2 = 0.3 and 
Pmaxl ' Po = 1-05. The structure of this mode will be compared in the following section with the results of numerical 
simulations. 



4. Numerical simulations 

In this section, we present the results of numerical simulations performed with ZEUS, NIRVANA and RAMSES. The 
goals are first to compare the properties of the instability during the linear phase of the parametric instability with 
the analytical theory presented above. Then we investigate the stability of strongly nonlinear waves. We also study the 
evolution of the instability during its nonlinear phase. In the next section we focus on waves for which (c/U) 2 = 0.3. 

4.1. Small amplitude waves 

We first study the case of a background wave with a small amplitude such that p m ax/po — 1-05. The properties of 
the unstable modes expected in that case were derived in section 13.41 From the solution of Eq. (|14[) we find that 
T w — 9.6H. The analytical results (sec figured]) show that the most unstable modes have the largest wavenumbers. 
However, because of their small scale structure, they will be difficult to represent numerically. In order to illustrate 
the properties of the parametric instability, we isolate the mode for which n = 4. This is achieved by adjusting the 
extent of the vertical domain appropriately. The results of section 13^41 indicate that X z = 2.77H and a — 0.06 for that 
mode. 

To follow its evolution, we performed 2D numerical simulations in the (x, z) plane with ZEUS, NIRVANA and 
RAMSES. For the n — 4 mode to develop, we chose a computational box that covers the range [— T w /2, T w /2] in x 
and [—X z /2, X z /2) in z (so that the radial extent matches the wavelength of the background wave while the vertical 
extent matches the vertical wavelength of the unstable mode). We used periodic boundary conditions in both x and z. 
Of course, other modes, compatible with the boundary conditions (i.e. those which have a ratio of box size to vertical 
wavelength that is a large integer) could in principle grow as well. However, very small scales modes will be strongly 
damped by numerical diffusion (see appendix^]). In practise, we found that they did not develop at the resolutions we 
investigated. We also found that an accurate description of the n = 4 mode requires roughly 64 cells per wavelength. 
This translates into a resolution (N X ,N Z ) — (320,64). For smaller resolutions, the growth rate is reduced, although 
the evolution is qualitatively similar. 

At time t = 0, the density and velocity components are set equal to the solutions of the system of differential 
equations given by Eq. (f2"2"]l and Eq. (|2"3")l . A perturbation is also added to each component of the velocity. If this 
perturbation was random, then a combination of modes with vertical wavenumbers +k z and —k z would grow with 
relative amplitudes depending on the initial conditions. This would make the comparison with the linear theory less 
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Fig. 6. The time history of the maximum value of the vertical velocity (left panel) and of the background wave 
amplitude (p m ax/po ~ 1) (right panel) as obtained with ZEUS (solid line), NIRVANA (dot-dashed line) and RAMSES 
(dashed line). The dotted line shown in the left panel indicates the behaviour expected from linear theory. 



illustrative. Rather, we chose to add a perturbation that is close to the mode with wavenumber +k z by adding a 
velocity perturbation of magnitude Sv such that 

Sv = Svq[— sm(nk x , w x — k z z) + sin((n + l)k x , w x — k z z)], (57) 

where we took <5t>o = 10~ 4 c. In practise, we found this procedure prevented the mode with wavenumber — k z from 
growing during the course of our simulations. 

The time history of the maximum value of the vertical velocity is plotted in the left panel of figure [5] The solid line 
shows the results obtained with ZEUS, while the dot-dashed and dashed line respectively correspond to the results 
obtained with NIRVANA and RAMSES. All curves indicate an exponential growth at early times. The evolution 
expected from linear theory is also plotted with the dotted line. There is very good agreement between this and the 
numerical results, which themselves are in very good agreement. Results from all three codes show that the linear 
instability saturates after about 100 wave periods, at which point the maximum vertical velocity has reached about 
10% of the sound speed. This is of the order of the radial velocity associated with the background wave. This indicates 
that the linear instability saturates when the amplitude of the velocities of the unstable inertial modes becomes large 
enough to perturb the structure of the background wave. During the nonlinear phase, the maximum value of v z at 
particular time is first found to oscillate before it saturates (as obtained with ZEUS) or slowly decays (as obtained with 
RAMSES). The oscillation is due to the linear instability recurring quasi periodically between t = 100 and t = 400, 
before the structure of the underlying wave is sufficiently modified. Finally, we comment that the different evolution 
of the simulations in the different codes is due to their different dissipative properties. In order to illustrate the effect 
of the instability on the wave structure, we plot the time history of its amplitude on the right hand side panel of figure 
[5] (with the same conventions as used on the left hand side panel). Once again, there is a good agreement between the 
codes during the linear phase of the evolution and during the recurring instability phase. RAMSES displays a larger 
amplitude than ZEUS and NIRVANA during that time. An additional run performed with NIRVANA using twice 
the resolution gave similar results to those of RAMSES, suggesting that this difference is due to the larger numerical 
dissipation in the finite difference schemes. At late times, all codes find that the wave amplitude is reduced to smaller 
and smaller amplitudes, although there are differences due to the different small scale dissipative properties. 

Figure [7] gives plots of the structure of the eigenmode obtained with the three different codes. The top panel shows 
the density distribution, normalised by the mean density, in the (r, z) plane after 50 wave periods, as obtained with 
ZEUS. Being in the linear phase, it is hardly modified compared to the initial profile. The next three panels, from top 
to bottom, represent the vertical velocity (normalised by the sound speed c) respectively with ZEUS, NIRVANA and 
RAMSES. They look almost identical to each other and extremely similar to the expected theoretical result shown in 
figure [5] 
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Fig. 7. Snapshots showing the structure of the parametric instability in the (r, z) plane in the case of a linear back- 
ground wave whose parameter are (c/U) 2 = 0.3 and p ma x/ Po — 1-05. The first two panels, both obtained with ZEUS, 
respectively show the density (normalised by the mean density po) and vertical velocity distribution (normalised by 
the sound speed). The third and forth panels show snapshots of the vertical velocity obtained with NIRVANA and 
RAMSES respectively. Both are very similar to the ZEUS results. All the plots are obtained at t = 50. 



Taken together, all the results presented here demonstrate a very good agreement between the numerical and 
analytical results. In the next section we investigate the properties of the parametric instability that destabilizes waves 
of larger amplitude. 



4.2. Larger amplitude waves 

Using ZEUS, we followed the evolution of waves of different amplitude. We used the same set up (resolution, initial 
perturbation) as for the small amplitude wave described above. In each case, the size of the computational box is 
tuned to match the background wave period in the radial direction and the wavelength of the n = 4 mode in the 
vertical direction. The results are illustrated in figure [8] We show the growth rates u of the instability measured during 
the linear phase for wave amplitudes such that p ma x/Po = 1-05, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7 and 1.8. They can be 
compared with the growth rates computed according to Eq. (f5"2l . To calculate these, the amplitude e that appears 
in Eq. (fB"2"|) is taken to be equal to the amplitude of the Fourier mode with wavelength equal to the radial width of 
the box. Figure [5] shows that there is very good agreement between the results of the numerical simulation and the 
analytical results for wave amplitudes smaller than p ma x/po = 1-5. For larger amplitudes, the analytical estimates 
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Fig. 8. Growth rate of the n = 4 mode of the parametric instability for different background wave amplitudes that 
are characterized by (c/U) 2 = 0.3. The stars are determined by fitting the results of numerical simulations obtained 
with ZEUS. The circles are obtained using the results of the linear analysis (using Eq. (|52|) . the parameter e is taken 
to be the amplitude of the Fourier mode whose wavelength equals that of the wave). 




Fig. 9. Same as figure [Jj but for a large amplitude background wave such that p m ax/pa = 1-8. The snapshots are 
obtained at time t = 10. Note the distorted structure of the eigenmode close to the density maximum. It results from 
the strong nonlinearities of the background wave at this location. 



tend to underestimate the growth rate slightly because of the strongly nonlinear nature of the wave. This nonlinearity 
modifies the structure of the eigenmode, as shown on figure [9] for the case p m ax/po — 1-8. This figure is similar to 
figure [7) the top panel shows the density distribution in the (r, z) plane and the bottom panel represents the vertical 
velocity. Both were obtained with ZEUS at time t = 10, well within the linear phase of the instability. A strong cusp, 
already noted in figure [2j is seen almost undisturbed on the upper snapshot. The lower one shows the structure on 
the eigenmode. It is similar to the plots of figure [JJ although its pattern is distorted close to the wave cusp. This is 
probably affecting the growth rate of the mode as indicated in figure [5] 

The long term evolution of the instability in this strongly nonlinear regime is further illustrated by figure ITDl which 
is the equivalent of figure [H ZEUS results are plotted using the solid line while RAMSES results are plotted with 
the dashed line. Both are almost identical and results obtained using NIRVANA displays a very similar behaviour. 
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Fig. 10. Same as figure but for a large amplitude background wave such that p m ax/po — 1-8 (only the ZEUS 
and RAMSES results are represented, respectively with the solid and dashed lines, while NIRVANA results are very 
similar). Note the faster growth of the instability compared to the smaller amplitude case. 
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Fig. 11. Structure of the flow during the nonlinear stage of the parametric instability. Left panels correspond to a 
background wave amplitude p max /po = 1-05 and (c/U) 2 = 0.3. The upper figure plots p/po in the (r, z) plane, while 
the lower figure shows v z /c. Both are plotted at t = 500. The right hand side panels shows the same quantities at 
t — 27 for a background wave amplitude p m ax/ Po=l-8 and (c/U) 2 = 0.3. 



The evolution of the instability is qualitatively similar to the small amplitude wave case, although the timescale for 
the instability to saturate is much shorter, being 20 wave periods when p m ax/po = 1-8 as compared to 100 wave 
periods when p m ax/po — 1-05 (once again, the linear instability saturates when the perturbed velocities reach that of 
the background wave, which happens to be close to the speed of sound in the large amplitude case). In addition, no 
recurring instability phase is found before the amplitude of the wave starts to decrease. This is better illustrated in 
figure QT] The left panels show the density structure (upper panels) and vertical velocity (lower panels) obtained in the 
small amplitude case (p max / Po=l-Q5) at time t = 500. Even at those late times, although smaller scale disturbances 
are involved, the regular pattern characteristic of the parametric instability is still visible. This is in contrast with the 
snapshots obtained in the large amplitude case (p m ax/po=l-8) at time t = 27 and shown in the right panels of figure 
ITT] The nature of the flow is far more disorganized and looks more like turbulence in that case. Finally, we want to 
comment that the system eventually goes back to equilibrium, with the wave amplitude being strongly damped. This 
might seem to disagree with figures [6] and [TU1 where the vertical velocity saturates and shows very little decrease. This 
is because a vertical flow of alternating direction sets in at large times due to the periodic boundary conditions in 
the vertical direction and damps only very slowly. v z being only a function of x, this flow simply superposes to the 
equilibrium without disturbing it. We emphasize here that this is only an artifact of the numerical setup. 
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5. Stratified case 

The results of the previous section give confidence that the nature of the instability is properly captured by the 
analytical work presented in section [3] However, this analysis neglected the disk density vertical stratification, which 
is the focus of the present section. 



5.1. Inertia I wave properties 

The properties of inertial waves in stratified isothermal disks have been studied bv lLubow fc Pringld (|l993h . We recall 



these below for the case 7 = 1 in their notation, as this is the case of interest here. Let v z (x,z) be the Eulerian 
perturbation of the velocity in the vertical direction. Because of vertical stratification, the z dependence in v z (x, z) is 
no longer harmonic as it was in the unstratified case. Instead, it can be shown to be of the form 

v z {x,y,z) = v z {z) exp[i(nk x , w x - u n t)} , (58) 

where v z (z) is a polynomial in z of the form: 

v~ z (z) = 2 c m z m . (59) 

m=0 

Here n z is given by 

where the frequency of inertial waves, which constitute the low frequency branch of this relation, is now denoted by 
u>n,n z since it depends on the two integ ers n and n z . In addition, the coefficients c m satisfy a recurrence relation, given 
by Eq. (49) of iLubow fc Pringle! (|l99sh . In the large n z limit, v z (z) oscillates with a period 



2ttH , . 

A z> „, ~ -— (61) 



in the vertical direction in the vicinity of the equatorial plane. Because of separability, two inertial waves with the 
same n z but different values of n have the same vertical structure. They can therefore interact together in the presence 
of the background nonlinear sound wave when they satisfy the same resonant condition as in the unstratified case: 



(62) 



As for the unstratified case, the spectral frequencies 0J n ,n z are dense in the interval [0,17] (this can be demonstrated 
by showing that the difference u> n +i,n z — u n ,n z tends to zero as n z tends to infinity): there will always be an infinite 
number of values of n z such that Eq. |62|) is satisfied to a specified accuracy. But note that the modes are in fact 
denser in the stratified case than in the unstratified case. This is because n z appears linearly in Eq. (|62p whereas the 
corresponding expression in the unstratified case would contain n\ (see Eq. (J5SJ). This has the consequence that the 
separation between neighbouring modes is less in the former case. A physical explanation for this is that this results 
from there being a continuous range of density scale heights in the stratified disk model compared to only one length 
scale in the unstratified case provided by the size of the computational box. This makes parametric instability more 
readily realised in the stratified case. 

In a numerical simulation, however, only those modes for which the wavelength is big enough (or n small enough) 
to be well resolved will have a non vanishing growth rate. This limits the number of realisable modes to those that 
satisfy the resonant condition to within a given accuracy and are such that n < n max . 



5.2. Numerical simulations 

The simulations we describe in this section are the stratified analogs of the simulations presented in section [4] (we only 
report results obtained with ZEUS here as we have shown good agreement between the different numerical methods 
above). The wave velocity is identical, (c/U) 2 = 0.3, and we consider only two wave amplitudes defined such that 
Pmax/Po = 1-05 and p ma x/Po = 1-8. In our standard runs, the resolution is kept the same in the x direction, N x = 320, 
as is the extent of the computational domain, x E [—T w /2, T w /2]. In the vertical direction, the computational domain 
is extended to cover 4 scale heights on both sides of the equatorial plane and the resolution is consequently increased to 
193 grid points. The initial density is set to be in hydrostatic equilibrium in the vertical direction (this gives a Gaussian 
vertical profile for the density) while the profile of the variables in the radial direction is unchanged compared to the 
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Fig. 12. Time history of the vertical velocity in the stratified model, for the small amplitude wave p m ax/po = 1-05 
(left panel) and for the large amplitude wave p m ax/po = 1-8 (right panel). In both cases, the dotted line represents 
the exponential fit obtained during the linear phase of the instability. The growth rates that result are respectively 
a = 0.038 and 0.32. 



unstratified case. Finally, we use periodic boundary conditions in the radial direction and outflow boundary conditions 
in the vertical directions. At t = 0, all components of the velocity are randomly perturbed with an amplitude equal 
to 0.01% of the sound speed. Note that this random initial condition does not break the symmetry of the problem at 
t = 0. Therefore, we expect modes with wavenumbers k z and — k z to develop at the same time. 

As in the unstratified simulations, the parametric instability is seen to develop in both runs. This is illustrated in 
figure 1121 which shows the time history of the maximum value of the vertical velocity in both cases (as before, time 
is measured in units of the wave period). The left panel corresponds to the small amplitude wave, Pmax/po = 1-05, 
and the right panel to the large amplitude wave p ma x/po = 1-8. Similarly to the unstratified case, the growth of the 
instability is faster in the latter case. The growth rates, measured during the linear phase and shown on figure 
using a dotted line, are respectively a = 0.038 and a = 0.32. Both are comparable (although slightly smaller) than 
the growth rates obtained in section R~2l and plotted in figure [5] We also checked that the amplitude of the wave starts 
to decrease as the instability enters the nonlinear regime and that the system goes back to hydrostatic equilibrium. 
The nonzero velocities seen at large times in figure rT2] are due to motions in the disk upper layers which are possibly 
long wavelengths acoustic modes excited by the parametric instability. Because of the low inertia of the gas, they take 
very long to damp. The structure of the flow during the linear phase is illustrated in both simulations by showing 
snapshots of different quantities in the (x, z) plane in figure [TBI The left hand side snapshots correspond to the small 
amplitude wave case and the right hand side snapshots to the large amplitude wave case (note that both plots only 
cover the domain [-2.5H, 2.5H] in the vertical direction, ie only a fraction of the actual computational box). On both 
sides, the upper panel shows the distribution of the density during the linear phase, respectively at times t = 150 and 
t = 17. As in figure [TOl the cusp in the density is apparent in the latter case. The middle panels in figure [TBI represent 
the vertical velocity at the same times as the upper panels. The stripped structure of the flow characteristic of the 
parametric instability is obvious in both simulations. In the small amplitude wave case, the structure of the unstable 
eigenmode is regular while it is distorted in the neighbourhood of the cusp in the large amplitude wave case. This 
is again similar to the results obtained in the unstratified simulations. Also similar is the outcome of the instability 
during the nonlinear regime. This is illustrated by the lower panels of figure [TBI which show the structure of the 
vertical velocity in both cases, respectively at times t = 450 and t = 45 (i.e. well into the nonlinear regime). The small 
amplitude wave simulation still displays a very regular structure while the entire flow broke down into turbulence in 
the large amplitude wave case. In conclusion, all of the results obtained in these stratified simulations are very similar 
to the unstratified results of section [TJ 

One may wonder, however, what determines the particular mode(s) that grow in these stratified simulations. Unlike 
the unstratified case, no vertical scale is determined by the size of the computational box. To get some insight into 
this question, we now compare the results of the small amplitude wave simulation to the linear analysis (performed 
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Fig. 13. Snapshots of the density and velocity in the (x, z) plane showing the structure of the parametric instability 
in a stratified disk. The panels on the left correspond to the small amplitude wave, p m ax/po — 1-05, while the right 
hand side panels refer to the large amplitude wave Pmax/po — 1-8. On both sides, the top panels show the density 
(normalised by the mean density) during the linear phase (respectively at t = 150 and t = 17), the middle panels 
show the vertical velocity (normalised by the sound speed) at the same times and the bottom panels show the vertical 
velocity during the nonlinear phase (respectively at t = 450 and t — 45). 




Fig. 14. Vertical velocity in the [x, z) plane resulting from the n = 9 mode only, calculated according to the linear 
results of section [3] in the case of the small amplitude wave p m ax/po = 1-05. It is seen to be in good agreement with 
the numerical simulations (see the snapshots given by figure I13[) . 




Fig. 15. Snapshot of the vertical velocity at time t — 200 in the low resolution run: (N x , N z ) — (160, 93) (left panel). 
The color scale spans the interval [—0.15c, 0.15c]. Note the larger scale of the unstable eigenmode that grows in that 
case as compared to the standard run with twice that resolution. It is compared on the right panel with the structure 
of the n = 4 eigenmode as calculated according to the results of section l3~4l (recall that both n and n + 1 contribute). 
Good agreement is found between the numerical and analytical results 
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Mode number 


n 


n z 


precision 


1 


9 


21 


2.3 x 1(T 3 


2 


10 


26 


4.6 x 10 -3 


3 


7 


13 


6.3 x 1CT 3 


4 


10 


25 


7.6 x 1CT 3 


5 


8 


17 


8.2 x 1(T 3 


6 


4 


4 


8.9 x 10" 3 



Table 1. Modes ranked according to the precision 9 (see its definition in the text) within which they satify the resonant 
condition given by Eq. (|62p . For each mode, labelled according to the first column, 9 is given in the last column, while 
the second and third columns respectively give the values of n and n z . The parameters of the model are such that 
(c/U) 2 = 0.3 and the background wave has an amplitude pj p ma x = 1-05. 



in the unstratified box) presented in section 13.41 We expect the latter to agree with the numerical results when large 
wavenumber modes are excited. Indeed, the instability is truly local in that case and should not be affected by the 
vertical density stratification. In figure [T3] the dominant mode that emerges in the simulations seems to be such that 
n ~ 8 to 10. As an example, the structure of the flow resulting from the mode n = 9 is illustrated on figure [14] 
(note that k z is determined as the value that leads to the satisfaction of the parametric resonance condition and both 
modes with wavenumber +k z and — k z are included). It is roughly in agreement with the left hand side snapshots of 
figure 1131 although some signs of the vertical stratification are seen in the latter through an increase of the velocity 
at large z. Figure IT51 also displays a signature of large scale variations in the radial direction, possibly indicating that 
a smaller n mode may be growing at the same time (see below). The linear analysis indicates that the n = 9 mode 
should have a wavelength A z = 1.34_ff, in rough agreement with the simulation. However, it predicts a growth rate 
vim = 0.062, larger than the numerically determined value a nurn — 0.038. This discrepancy is due partly to numerical 
damping reducing the growth r ate. It is also du e to the resonant condition being actually only marginally satisfied, 
which reduces the growth rate (lGoodmanlll993h . We quantify that accuracy using the parameter 9 which we define 
through 

q ^n,n z T" bJn+l,nz ^ (63) 

In the notation of section 1331 9 = Su/uj. Eq. ([53"]) states that the growth rate is unchanged when 9 <C <r/7r, while 
the instability is suppressed when 9 > <j/tt. According to the linear results of section l3"^il a/ir ~ 0.01 for the small 
amplitude wave parameters we are studying here. We therefore expect that only those modes that have 9 < 0.01 will 
grow with a significant growth rate. Using Eq. (f60|) . we searched for all the modes having n < 10 and 9 < 0.01. The 
modes satisfying these criteria are listed in Table [TJ The mode with the smallest value of 9 has (n, n z ) = (9, 21) and 
probably corresponds to the mode described above and seen in the simulations. In fact, Eq. (j6Tj) gives A 2i „ 2 = 1.37 H 
for that mode, which is in good agreement both with the results of the linear analysis mentioned above and to the 
results of the numerical simulations. All the other modes, with n £ [7, 10] and increasing values of 9 probably contribute 
in a complicated way to the growth of the parametric instability. Table 1 also reports the existence of a mode with 
(n,n z ) = (4,4) and 9 = 8.9 x 10 -3 which probably accounts for the large scale structure seen in figure [TBI As the 
resolution increases, we expect more modes with smaller wavelength to become important and to contribute to the 
structure of the flow. This is also expected as a consequence of numerical dissipation acting on the smallest scales 
(see appendix [A]) . To illustrate these points, we performed a simulation with the resolution decreased by a factor of 
two: (N X ,N Z ) — (160,93). Everything else is kept identical. The structure of the flow in that case is illustrated on 
the left panel of figure [15] This is again a snapshot of the vertical velocity in the (x, z) plane at time t = 200. In 
agreement with the previous discussion, the growth of the parametric instability is dominated by an eigenmode whose 
wavelength is larger than the high resolution run. In fact, its structure probably results from the growth of the n = 4 
mode listed in table [TJ To check that hypothesis, the flow is compared in the right panel with the structure of the 
n = 4 mode as computed according to the unstratified linear analysis. The agreement between the analytical results 
and the numerical simulations is very good. 



6. Conclusion 

In this paper, we have derived and solved a simple ODE which describes the structure of nonlinear density waves 
propagating in accretion disks. Such an analysis is important for understanding the properties, stability and dissipation 
of waves that propagate in disks in many situations: tidal excitation by a companion, waves generated by the disk 
self-gravity when it is massive enough or even waves generated by MHD turbulence resulting from the MRI. For 
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the analysis to remain tractable, we adopted in the present paper the shearing box approximation and restricted the 
equation of state to be isothermal and the wave to have no variation in the direction of the mean shear flow ie. to 
be axisymmetric. Because the waves were free and stationary in a frame moving with the phase velocity, shocks that 
would require a continual forcing or energy input were not considered here. 

Using the solution of the ODE that describes the propagation of waves of arbitrary amplitudes to provide ini- 
tial conditions, we followed their subsequent evolution in ID using different numerical methods. A Godunov code, 
RAMSES, was found to give very accurate results. Finite difference codes such as ZEUS or NIRVANA were found to 
diverge from the correct solution or develop small scale oscillations on long timescales. These problems were cured in 
our case by decreasing the Courant number to C — 0.1. However, we comment that reducing the timestep only delays 
the onset of the problems which were always found to develop provided the waves were evolved for long enough. 

We also studied the stability of the waves. When their frequency is in the range [fl, 20], they can interact resonantly 
with two linear inertial waves and undergo a parametric instability. In this context we note that th is frequency range 
may be increased in stratified disks if 7 7^ 1 for perturbations allowing the existence of g modes ( Lubow fc Pringlel 



19931) . A wider class of waves would thus become candidates for parametric instability. Neglecting vertical density 
stratification, we derived an expression for the growth rate when the amplitude of the background wave is small. This 
was found to be in good agreement with 2D numerical simulations using RAMSES and ZEUS or NIRVANA (provided 
that in the latter cases, the timestep was tuned as described above). Both unstratified and stratified simulations 
were carried out. We found that the precise structure of the mode that grows depends on the resolution we used. 
This is because we used no physical dissipation and rather relied on the resolution dependent numerical viscosity to 
damp small scale modes. Using these simulations, we could study the nonlinear evolution of the instability. For small 
amplitude waves the motion remained ordered but large amplitude waves broke down to a turbulent state as might 
be expected from consideration of the related elliptical instability in rotating flows. In all cases, we found that the 
amplitude of the wave was damped as the system relaxed back to a stationary state. Of course, this is because no 
external energy that could keep the instability going is provided. This differs from realistic situations that can occur, 
such as when the disk is tidally forced by a companion. Then the background wave will be constantly forced. In such 
cases, the external forcing will keep driving the instability which we expect to result in a complicated quasi-steady 
state. We will study forced waves, which also allow the incorporation of shocks, in a future paper. 
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Appendix A: The effect of a physical dissipation 

Here we consider the effect of a small constant kinematic viscosity v on the linear growth of the parametric instability. 
To do so, we study the properties of unperturbed inertial waves (see section 13.21) when v is nonzero. Note however 
that we neglect the effect of the dissipation on the background wave because its scale is much larger. In the large 
wavenumber limit, we may neglect the effect of possible stratification, the inertial modes are nearly incompressible 
and Eq. (|36p and ([37| can be written as 



— £,4- Sl 2 £„ = c 2 — 4 c 2 — — 4- iA7 2 U— £, 



dt dx ) ^ z dz 2 dxdz \ dt dx ) ^ x ' ^ ^ 

when viscosity is taken into account. Using Eq. (|38[) . these two equations lead to: 

[<4 - n 2 k 2 x w c 2 - ivk 2 u n ]Q + nk x , w k z c 2 C = , (A.3) 

nk x , w k z c 2 Q + K 2 - k 2 z c 2 - ivk 2 uj n \& = , (A.4) 

where uj n — U)i >n — nu> (note that we assumed that k 2 c 2 3> O 2 in writing the above equations). Neglecting terms in 
order v 2 , these equations are in fact identical to their inviscid counterparts, Eq. (|3"9")) and (f4U)l. provided w„ is replaced 
by u>n,visc = <^n — ivk 2 /2. This indicates that inertial waves are damped on a timescale 

Tdamp — ■ (A-"-*) 

vk ,L 

There will be a corresponding reduction in the growth rate of the parametric instability. In numerical simulations 
performed with Eulerian codes, we assume that the numerical dissipation has an effect similar to that describe above, 
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with a "numerical" viscosity v num . For a given resolution, v num is fixed and Eq. (|A.5[) shows that modes with larger 
wavenumbers are damped more rapidly, in agreement with our findings. 

Another application of this result is to estimate a numerical Reynolds number Re for our simulations. Re is given 

by 

cH 

Re = . (A.6) 

"num 

In section 14.11 we found that the linear growth rate of the n = 4 mode is unaffected by numerical dissipation, which 
indicates that 

Tdamp 3> ■ (A. 7) 

a 

This gives an upper limit for the numerical viscosity of the codes and a lower limit for the Reynolds number of our 
simulations. Introduc ing the parameters of the simulation in Eq. (|A.7[) . we found Re > 10 3 . This is consistent with 
the value reported bv lBalbus et al. I (|l996h for the Reynolds number of Eulerian codes with similar resolutions. 
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